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IX 


1.  Summary 


Explosions  are  often  conducted  in  complexes  with  chambers,  tunnels,  and  shafts 
used  for  access  and  instrumentation.  These  structures  and  sharp  topographic  features  can 
act  as  strong  scatterers  of  seismic  waves  and  complicate  the  radiation  patterns  from 
explosions.  In  extreme  cases  they  could  affect  discrimination  between  earthquakes  and 
explosions.  The  objectives  of  this  research  are  (1)  to  study  the  effects  of  these  near¬ 
source  scatterers  on  seismic  waves  radiated  from  explosions,  and  (2)  to  use  a  wavelet 
domain  based  moment-tensor  inversion  scheme  to  determine  “explosive”  and  “multi¬ 
couple”  components  of  the  source. 

To  calculate  seismograms  from  an  explosion  near  strong  scatterers  we  use  a  new 
finite  difference  algorithm  that  accomplishes  variable  griding  by  coordinate 
transformation  or  “stretching.”  This  method  provides  excellent  numerical  stability  while 
increasing  computational  efficiency.  It  is  capable  of  3-D  calculations  for  sources  near 
strong  scatterers  in  heterogeneous  media.  Seismograms  are  calculated  to  determine 
effects  of  various  scatterers  on  seismic  radiation  patterns.  The  code  is  developed  for 
realistic  earth  models  including  (1)  free  surface,  (2)  layered  structure,  (3)  surface 
topography,  and  (4)  seismic  attenuation.  In  addition,  a  perfectly  matched  layer  (PML) 
was  incorporated  into  the  finite-difference  code  to  improve  the  absorption  at  the 
boundaries  and  for  saving  memory.  Forward  modelings  using  the  3-D  finite-difference 
code  are  conducted  with  an  explosive  source  and  tunnel  in  a  full  space  and  a  layered  half¬ 
space.  Calculations  are  carried  out  for  each  case:  (1)  a  reference  model  without  a  tunnel, 
and  (2)  a  finite  length  horizontal  tunnel  included.  The  calculations  show  P  to  P  and  P  to 
S  scattering  and  a  complicated  radiation  pattern  of  the  wave  field.  P  to  S  scattering  is 
strong  and  the  tunnel  acts  as  a  virtual  shear  wave  source.  For  the  half-space  model, 
surface  waves  dominate  the  seismograms  and  significant  SH  waves  are  generated  by  the 
presence  of  the  tunnel  because  of  shallow  depth.  We  show  the  capacity  of  our  code  in 
modeling  wave  scattering  due  to  various  topographical  features.  The  sharp  comers  of  a 
mesa  act  as  strong  scatters  of  seismic  waves,  and  a  smooth  hill  scatters  the  waves  much 
less. 
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We  use  synthetic  seismograms  to  test  the  performance  of  moment  tensor 
inversion  and  its  ability  to  separate  the  volumetric  and  shear  components  of  the  source. 
With  good  azimuthal  coverage,  the  moment  tensor  shows  significant  shear  components  in 
the  presence  of  a  scatterer  when  surface  waves  dominate  the  seismogram.  The  use  of 
wavelets  for  moment-tensor  inversion  has  advantages  because  it  can  use  all  or  part  of  a 
seismogram,  and  has  effective  de-noising  capability.  The  method  can  be  used  to  analyze 
data  from  earthquakes  and  explosions  for  determining  isotropic  (explosion)  and  multi¬ 
couple  (earthquake)  components  of  the  source.  This  is  valuable  for  seismic 
discrimination. 

We  test  the  applicability  of  Time  Reversed  Acoustics  (TRA)  approach  for 
detecting  the  presence  of  a  tunnel  near  the  source.  In  this  approach,  the  recorded 
seismograms  are  time-reversed  and  sent  back  into  the  earth  at  each  station.  The  back- 
propagating  wavefields  focus  at  the  source.  The  P  wave  focuses  strongly  at  the  explosion 
while  the  S  wave  at  the  tunnel.  TRA  has  great  potential  for  determining  the  seismic 
source  properties. 
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2.  Mechanisms  on  Scattering  due  to  an  Explosive  Source 


2.1  Observation  on  Transverse  wave  generations 

An  explosive  source  in  a  laterally  homogeneous  layered  half-space  generates  P, 
SV,  and  Rayleigh  waves,  but  not  SH  or  Love  waves.  However,  seismograms  from  a 
large  number  of  explosions  show  nonisotropic  radiation  patterns  for  P  and  Rayleigh 
waves,  and  prominent  SH  and  Love  waves.  Various  mechanisms  have  been  proposed  to 
explain  the  generation  of  these  transverse  waves. 

2.2  Objectives  of  This  Research 

In  this  project,  we  will  study  the  role  of  scattering  from  near  source  structures 
such  as  tunnels,  shafts,  adits,  and  surface  topography  for  contributing  to  the  generation  of 
complicated  radiation  patterns  and  SH  waves  from  explosions.  In  addition,  we  will 
utilize  a  wavelet  domain  moment  tensor  inversion  to  separate  the  isotropic  (i.e., 
explosion)  and  multi-couple  components  of  the  source  as  an  aid  to  seismic 
discrimination. 

Medium  properties  near  the  source,  heterogeneities  along  the  propagation  paths, 
and  scattering  near  the  receivers  (Johnson,  1995)  all  contribute  to  complexities  of  seismic 
waves  recorded  from  explosions.  Path  and  near  receiver  effects  are  unique  to  each 
receiver.  The  near-source  complexities,  however,  affect  all  seismograms  and  lead  to 
complicated  source  functions  that  affect  the  important  task  of  discriminating  between 
earthquakes  and  explosions. 

The  near-source  contribution  to  nonisotropic  radiation  of  P,  Rayleigh  waves  and 
SH  wave  generation  by  an  explosion  can  be  attributed  to  three  mechanisms.  First,  is  the 
tectonic  strain  energy  released  by  the  explosion,  resulting  in  a  composite  source 
consisting  of  an  explosion  and  a  double  couple  (e.g.,  Archambeau  and  Wqmmiw,  1970; 
Toksbz  et  al.,  1965;  Toksoz  and  Kehrer,  1972;  Wallace  el  al .,  1983,  1985;  Priestley  et  al., 
1990;  Schlittenhardt,  1991;  Stump  et  al.,  1994;  Patton  and  Taylor,  1995).  Tectonic  strain 
release,  widely  observed  for  explosions  in  competent  rock,  has  a  greater  contribution  at 
low  frequencies  than  at  high  frequencies. 
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The  second  factor  that  contributes  to  nonisotropic  seismic  radiation  pattern  from 
explosive  sources  is  the  shape  of  the  explosion  cavity,  and  the  location  of  explosive 
source  in  that  cavity  (Rial,  and  Moran,  1986;  Stevens  et  al.,  1991;  Zhao  and  Harkrider, 
1992;  Ben-Menahem  and  Mikhailov,  1995;  Gibson  et  al.,  1996;  Ben-Menahem,  1997; 
Imhof  and  Toksoz,  2002;).  The  contribution  of  this  mechanism  to  nonisotropic  radiation 
depends  on  the  aspect  ratio  of  the  cavity  (i.e.,  ratio  of  length  to  diameter),  and  placement 
of  the  source  within  the  cavity  (Gibson  et  al.,  1996).  To  have  a  significant  effect,  a  very 
large  aspect  ratio  or  a  source  placed  in  one  end  of  a  long  cylindrical  cavity  was  required. 
Stevens  et  al.  (1991)  examined  the  effects  of  an  ellipsoidal  cavity.  Their  conclusion  was 
that  at  frequencies  of  1  Hz  or  less,  the  P-wave  radiation  pattern  was  nearly  isotropic  and 
the  S-wave  amplitudes  were  small. 

The  third  mechanism  contributing  to  complex  seismic  radiation  from  an  explosion 
is  the  near  source  scattering  (Gupta  et  al.,  1990;  Johnson,  1997;  Ben-Menahem,  1997; 
Imhof  and  Toksoz,  2002,  002).  This  mechanism,  not  studied  as  extensively,  could  be 
significant  when  strong  scatterers  are  present  near  the  source.  The  first  task  of  the 
proposed  research  is  this  near  source  scattering  mechanism. 

2.3.  Effect  of  Near  Source  Scattering  on  Seismic  Waves  from  Explosions 

Explosions  often  take  place  in  complexes  with  chambers,  shafts,  and  tunnels  used 
for  access  and  instrumentation  (Denny  and  Stull,  1994;  Muiphy  et  al.,  1994).  Cavities 
and  tunnels  can  act  as  strong  scatterers  and  efficiently  convert  different  wave-types  into 
each  other  (Tadeu  et  al.,  1996).  The  scatterers  act  as  secondary  sources  with  difference 
radiation  patterns  and  different  distance  dependencies  than  the  primary  waves.  A  large 
amount  of  S  wave  energy  can  be  generated  by  this  mechanism  by  an  explosion  (Leavy, 
1993;  Ben-Menahem,  1997;  Imhof  and  Toksoz,  2000,  2002). 

Imhof  and  Toksoz  (2000,  2002)  modeled  the  effects  of  scattering  from 
cylindrical  cavities  near  a  shot  using  a  multiple-multipole  approximation.  Their  results 
showed  scattering  affected  the  P-wave  radiation  pattern  and  generated  significant  S 
waves.  Because  of  the  limitation  of  the  method,  these  studies  were  limited  to  2D 
geometry,  with  infinitely  long  cylinders  for  sources  and  scatterers. 
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To  evaluate  the  effects  of  scattering  realistically,  one  needs  to  solve  the  3D 
problem  consisting  of  explosive  sources,  finite  dimension  scatterers  of  various  shapes 
(e.g.,  cylinders),  and  a  free  surface.  This  is  a  formidable  numerical  problem,  but  can  be 
done  with  an  appropriate  numerical  approach.  For  this  we  will  use  a  new  finite 
difference  algorithm  described  below. 

3.  Application  of  Finite  Difference  Modeling  in  Studying  Wave 
Propagation 


3.1  Finite  Difference  with  Coordinate  Stretching 

The  finite  difference  time  domain  method  (FDTD)  is  one  of  the  most  widely  used 
tools  to  simulate  wave  propagations  in  2-D  and  3-D  elastic  media  with  general  spatial 
variations  of  elastic  properties.  While  efficiency  can  be  achieved  by  sampling  the 
physical  space  adaptively  with  a  variable  grid,  this  benefit  may  be  offset  by  problems 
introduced  by  change  in  grid  size.  Wave  distortion  or  numerical  reflections  may  occur 
due  to  phase  change  at  the  interface  of  two  grids  (Browning  et  al.,  1973).  Numerical 
reflections  can  be  eliminated  by  taking  extra  averages  of  solutions  at  the  grid  interface, 
however  wave  distortion  cannot  always  be  avoided.  It  is  also  hard  to  match  and 
implement  finite  difference  approximations  at  the  interface  of  different  grids  (Lilia,  1997; 
Hayashi,  1999).  Introducing  a  change  in  grid  spacing  may  adversely  affect  the  formal 
truncation  error  and  the  stability  of  the  system  (Crowder  and  Dalton,  1971).  By  counting 
phase  shift  on  a  nonuniform  spacing  grid,  a  3D  elastic  FDTD  approach  was  developed  to 
avoid  numerical  reflections  and  reduce  wave  distortions  (Pitarka,  1999). 

Pitarka’s  approach  requires  solving  a  linear  system  before  conducting  the  FDTD 
calculation.  The  linear  system  must  be  solved  again  when  changing  the  difference 
operator  (e.g.,  from  4th  to  8th  order).  To  avoid  doing  this,  we  apply  a  coordinate 
transformation,  or  stretching,  to  achieve  effective  variable  griding  in  the  physical  domain 
and  uniform  griding  in  the  transformed  domain.  Because  the  computation  is  carried  out 
on  a  uniform  grid,  the  formal  truncation  error  and  the  stability  properties  of  the  finite 
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difference  computation  are  preserved.  In  addition,  numerical  implementation  becomes 
easier  and  more  flexible. 

A  further  increase  in  computational  efficiency  and  accuracy  is  achieved  by 
employing  a  wavelet-based  difference  operator.  At  the  same  accuracy,  the  wavelet-based 
FDTD  scheme  allows  eight  times  more  coarse  sampling  for  3-D  models  than  the  popular 
4th  order  method,  and  is  comparable  to  higher  order  (6th/8th)  schemes.  At 
discontinuities,  the  wavelet-based  approach  outperforms  the  2nd,  4th,  6th  and  8th  order 
Taylor’s  expansion-based  schemes,  especially  when  combined  with  variable  griding  in 
the  neighborhood  of  a  discontinuity. 

3.2  Test  of  Scattering  Due  to  a  Near-Source  Cavity  in  2-Dimension 

As  a  test  and  also  the  first  step  of  this  research,  strong  scattering  from  the  near 
source  cavity  is  simulated  with  the  FDTD  algorithm  described  above.  The  model  is  2400 
m  X  2400  m  in  the  x  and  y  direction,  and  infinite  in  the  z  direction.  An  explosion  source 
with  center  frequency  of  10  Hz  is  placed  in  the  center  of  the  model.  The  cavity  is  an 
infinite  cylinder  full  of  air,  and  has  a  radius  of  10  m.  The  distance  between  the  source 
and  the  cavity  is  100  m.  The  compressional  and  shear  velocities  are  2000  m/s  and  1000 
m/s,  respectively.  Grid  size  for  the  solid  domain  is  10  m  and  l  m  for  the  cavity  area.  A 
smooth  transition  zone  of  10  m  is  applied  between  the  fine  and  coarse  grid  area  to  further 
reduce  grid-induced  reflection.  The  size  of  the  constructed  mesh  is  284  x  284.  Figures 
la  and  lb  show  the  model  setup  and  griding  details  close  to  the  cavity. 

Two  snapshots  of  the  stress  component  (Txx)  are  shown  in  Figures  2a, b.  These 
show  the  interaction  of  the  waves  with  the  cylindrical  cavity.  The  cavity  acts  as  a 
secondary  source  (Figure  2b)  radiating  both  P  and  S  waves. 

Figures  3a, b  and  4a,b  show  seismograms  as  a  function  of  azimuth  at  two 
horizontal  planes  located  100  m  and  500  m  above  the  source,  respectively.  The  effect  of 
the  cavity  alters  the  azimuthal  isotropy  of  P-waves  radiation.  Also,  there  is  a  small 
azimuth  dependent  tangential  component  to  P-waves.  The  SH-wave  radiation  pattern 
(Figures  3b  and  4b)  is  more  complicated  than  that  of  a  double  couple. 

As  part  of  the  work,  we  upgraded  the  program  from  2.5  dimensions  to  three 
dimensions.  With  3D  calculations,  scatterers  with  different  shapes  and  orientations, 
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velocity  variation  in  the  medium,  and  free  surface  effects  can  be  modeled.  The  original 
program  was  only  capable  for  2.5  dimensions  because  of  limitation  of  computer 
resources.  A  new  computer  cluster  with  24  CPUs  facilitated  the  3D  computation. 
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Figure  la.  Model  setup  and  station  map. 
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Figure  lb.  Griding  details  close  to  the  cavity. 
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(a)  snapshot  at  0.37005  second 


(b)  snapshot  at  0.82005  second 

Figure  2.  Snapshot  at  (a)  an  early  time  when  the  cavity  start  interfere  the  source 
radiation  pattern  and  another  one  at  (b)  a  later  time  when  both  the  direct  and 
scattered  field  are  developed.  The  slice  is  in  the  source  plane. 


8 


Time  (s)  Time  (s) 


Time  (s) 

(a)  radial  particle  velocity 


-i  o  1 


Time  (ms) 

(b)  tangential  particle  velocity 


Figure  3.  Radial  and  tangential  components  of  particle  velocities  at  each  azimuthal 
direction.  Receivers  are  100  m  above  the  source  plane. 
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Figure  4.  Radial  and  tangential  components  of  particle  velocities  at  each  azimuthal 
direction.  Receivers  are  500  m  above  the  source  plane. 
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3.3  Scattering  Due  to  a  Near-Source  Tunnel  in  3-Dimension 

Strong  scattering  from  a  near-source  tunnel  is  simulated  with  the  FDTD  algorithm 
described  above.  Figure  5  shows  the  3-D  geometry  of  the  simulation  model  in  the 
presence  of  a  cylindrical  tunnel. 


z 


Figure  5.  The  geometry  diagram  of  the  model.  The  explosion  source  is  set  in  the 
middle  of  the  model,  coincident  with  the  origin  of  the  coordinates.  A 
cylindrical  tunnel  is  set  30  m  away  from  the  source,  with  radius  of  15  m  and 
length  of  50  m.  Its  symmetric  axis  is  parallel  to  the  y  axis.  The  receivers  are 
set  in  the  plane  200  m  above  the  source.  The  distance  between  the  neighbor 
receivers  is  100  m. 

Simulation  is  also  conducted  when  the  tunnel  is  absent.  The  infinite  formation  is 
modeled  by  a  finite  volume  of  3000  m  by  3000  m  by  1 100  in  the  x,  y,  and  z  directions. 
Absorbing  boundaries  are  placed  at  all  six  sides  of  the  model  to  eliminate  wave 
reflection.  An  explosion  source  with  a  center  frequency  of  40  Hz  is  excited  in  the  center 
of  the  model.  The  length  and  radius  of  the  cylindrical  tunnel  are  50  m  and  15  m, 
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respectively.  The  axis  of  the  tunnel  is  parallel  to  the  y  axis  of  the  Cartesian  coordinate. 
The  distance  between  the  source  and  the  tunnel  is  30  m.  The  formation  is  a  perfect 
Poisson  medium  with  compressional  and  shear  wave  velocities  and  density  of  4000  m/s, 
2300  m/s,  and  2500  kg/m3,  respectively.  Inside  the  tunnel,  we  assume  a  velocity  of  340 
m/s  and  a  compressed  air  density  of  200  kg/m3.  Grid  size  increases  from  1  m  near  the 
source  and  cylinder  tunnel  to  10  m  away  from  the  source  region.  Examples  of  calculated 
seismograms  are  shown  in  Figure  6.  The  receiver  plane  is  in  the  XY  plane  and  200  m 
above  the  source.  The  seismograms  are  along  a  line  parallel  to  the  X  axis.  Computations 
are  made  for  an  explosion  alone  (without  a  tunnel)  and  an  explosion  with  a  tunnel 
present. 

The  first  row  of  Figure  6  shows  the  x  and  z  components  of  the  motion  (velocity 
fields)  when  the  tunnel  is  absent.  As  expected,  only  P  waves  are  generated.  When  the 
tunnel  is  present  near  the  explosion  source,  strong  scattering  occurs.  The  seismograms  in 
the  second  row  of  Figure  6  show  that  strong  shear  waves  are  scattered  from  the  tunnel. 
To  demonstrate  the  effect  of  scattering,  we  subtract  the  seismograms  with  an  explosion 
only  from  those  with  the  explosion  plus  the  tunnel.  The  third  row  of  Figure  6  shows 
strong  scattered  waves.  These  are  both  P  and  S,  and  scattered  S  waves  are  larger  than  P 
waves.  SH  waves  are  also  observed,  but  their  amplitudes  are  much  smaller  than  the  x  and 
z  components. 

Divergence  and  curl  of  the  wavefield  separate  the  compressional  and  shear 
components  of  the  wavefield.  Figure  7  shows  the  snapshots  of  the  divergence  and  curl  of 
the  velocity  field  at  60  ms  on  the  XZ,  YZ  and  XY  planes  passing  through  the  origin  of 
the  coordinate  system.  One  can  clearly  see  the  complex  pattern  of  the  scattering.  The 
curl  is  zero  inside  the  tunnel.  The  snapshots  of  the  XZ  plane  and  XY  plane  reveal  the 
vertical  and  horizontal  cross-section  of  the  tunnel.  The  snapshot  along  the  YZ  plane, 
parallel  to  the  tunnel,  shows  the  complexity  introduced  by  the  finite  length  of  the  tunnel. 

Understanding  the  patterns  of  scattered  waves  will  help  us  choose  appropriate 
models  for  inversion  of  the  source  mechanism.  In  the  XZ  plane,  the  shear  wave  energy  is 
divided  into  four  quadrants,  and  they  are  in-phase  in  the  diagonal  quadrants.  Shear  waves 
in  the  half-space  above  and  below  the  XY  plane  arc  180°  out  of  phase.  Disturbance  of 
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the  P  wave  field  due  to  scattering  can  also  be  seen  in  the  divergence  field,  but  the 
scattering  patterns  are  not  as  clear  as  those  of  the  shear  waves. 


scattered  vx  scattered  vz 


Figure  6:  Comparison  of  the  x  and  z  components  of  the  velocity  fields  at  the  receiver 
plane  in  the  absence  and  the  presence  of  the  tunnel.  Respective  scattered 
velocity  fields  of  the  x  and  z  components  are  also  shown  in  the  last  row  of  this 
Figure. 


13 


-200 


-TOO 


200 


with  cylinder,  curl  of  velocity  field  at  60  ms 


YZ 


with  cylinder,  div  of  velocity  field  at  60  ms 


-200 


-100 


100 


2U0 


-200 


200 


distance  (m) 


200 


-100  0  100 
distance  (m) 


200 


Figure  7:  Snapshots  of  the  divergence  and  curl  of  the  velocity  field  from  an 
explosion  in  the  presence  of  a  near-source  tunnel.  The  snapshots  were  taken  at  60 
ms  after  the  explosion. 
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The  presence  of  several  tunnels  creates  very  complicated  scattered  waves.  Figure 
8  shows  one  example  of  such  a  model,  where  there  are  six  intersecting  tunnels,  30  meters 
below  the  surface.  Two  components  of  the  seismic  motion  are  shown. 
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Figure  8:  Snapshots  of  the  Vx  and  Vz  components  of  the  velocity  field  from 
an  explosion  in  the  presence  of  multiple  near-source  tunnels.  The  snapshots  were 
taken  at  24,95  ms  after  the  explosion. 


100 


15 


3.4  Scattering  Modeling  Considering  Free  Surface  and  Attenuation 

We  enhanced  our  modeling  ability  to  realistic  earth  models  that  include  (1)  planar 
free  surface,  (2)  layered  structure,  and  (3)  seismic  attenuation.  We  also  incorporated 
perfectly  matched  layer  (PML)  into  the  finite-difference  code  to  improve  the  absorption 
at  the  boundaries  and  for  memory  saving.  The  program  can  easily  handle  arbitrary 
topography  by  using  a  rotated-staggered  grid  scheme. 

3.4.1  Model  Description 

Strong  scattering  from  a  near-source  tunnel  is  simulated  with  the  FDTD  algorithm 
described  above.  Figure  9  shows  the  3-D  geometry  of  the  two-layer  model  with  a 
cylindrical  tunnel  in  the  first  layer.  Simulation  is  also  conducted  without  the  tunnel.  The 
dimensions  of  the  model  in  the  x,  y,  and  z  directions  are  6,000  m,  6,000  m,  and  400  m, 
respectively.  The  thicknesses  of  the  first  and  second  layers  are  300  and  100  meters, 
respectively.  PML  absorbing  boundaries  are  placed  at  five  sides  of  the  model  to  eliminate 
wave  reflection.  The  top  of  the  model  is  the  free  surface.  An  explosion  source  with  a 
center  frequency  of  10  Hz  is  placed  at  100  m  depth.  The  length  and  radius  of  the 
cylindrical  tunnel  are  50  m  and  15  m,  respectively.  The  axis  of  the  tunnel  is  parallel  to 
the  y  axis  of  the  Cartesian  coordinates.  The  distance  between  the  source  and  the  tunnel 
axis  is  55  m. 

The  elastic  properties  of  the  formations  and  tunnel  are  as  follows.  In  the  first 
layer,  the  formation  compressional  and  shear  wave  velocities  and  density  are  3000  m/s, 
1700  m/s,  and  2300  kg/m3,  respectively.  In  the  second  layer,  the  formation 
compressional  and  shear  wave  velocities  and  density  are  4000  m/s,  2300  m/s,  and  2800 
kg/m3,  respectively.  Inside  the  tunnel,  we  assume  an  acoustic  velocity  of  340  m/s  and  a 
compressed  air  density  of  200  kg/m3. 

3.4.2  Modeling  results — Synthetic  seismograms 

The  receiver  array  with  100  m  spacing  forms  a  dense  grid  on  the  free  surface.  The 
first  row  of  Figure  10  shows  the  x  and  z  components  of  the  seismograms  along  the  x  axis 
(velocity  fields)  when  the  tunnel  is  absent.  Strong  surface  waves  are  observed  in  both  the 
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x  and  z  components.  Due  to  multiple  reflections  between  the  free  surface  and  the  layer 
boundary,  P,  S,  and  surface  waves  are  generated  and  the  wavefield  becomes  extremely 
complicated.  When  the  tunnel  is  present  near  the  explosion  source,  strong  scattering 
occurs.  The  seismograms  in  the  second  row  of  Figure  10  show  that  strong  surface  waves 
are  scattered  from  the  tunnel.  To  demonstrate  the  effect  of  scattering,  we  subtract  the 
seismograms  with  an  explosion  only  from  those  with  the  explosion  plus  the  tunnel.  The 
third  row  of  Figure  10  shows  strong  scattered  waves.  These  are  surface,  P,  and  S  wave; 
scattered  surface  and  S  waves  are  much  larger  than  P  waves.  SH  waves  are  also  observed 
in  the  transverse  components. 

Figure  11a  shows  the  azimuthal  traces  when  a  tunnel  is  present  near  the 
explosion.  Figure  lib  shows  only  the  scattered  wavefields.  The  forward  and  backward 
scatterings  are  much  stronger  than  in  the  direction  parallel  to  the  tunnel.  However,  SH 
components  can  be  observed  at  all  azimuths. 
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Figure  9:  The  geometry  diagram  of  the  tunnel  model.  The  explosion  source  is 
located  at  100  meter  deep.  A  cylindrical  tunnel  is  set  30  m  away  from  the  source, 
with  radius  of  15  m  and  length  of  50  m.  Its  symmetric  axis  is  parallel  to  the  y  axis. 
The  receivers  are  set  on  the  free  surface.  The  distance  between  the  adjacent 
receivers  in  x  or  y  axis  direction  is  100  m. 
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Figure  10:  Comparison  of  the  x  and  z  components  of  the  velocity  fields  at  the 
free  surface  in  the  absence  and  the  presence  of  the  tunnel.  Respective  scattered 
velocity  fields  of  the  x  and  z  components  are  also  shown  in  the  last  row  of  this 
Figure. 
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Figure  11:  Three  components  of  the  recorded  wavefields  at  1km  offset  of  different 
azimuths,  a)  Free  surface  with  tunnel  near  the  explosion,  b)  Scattered  wave  only. 


3.5  A  rotated-staggered  grid  for  finite-difference  modeling 

We  implement  a  rotated-staggered  grid  (RSG)  for  finite-difference  modeling. 
Now  the  program  can  easily  handle  arbitrary  topography  by  using  a  rotated-staggered 
grid  scheme.  The  program  can  also  be  run  on  the  large  scale  Alliance  for  Computational 
Earth  Sciences  cluster.  It  also  adopts  a  perfectly  matched  layer  boundary  condition  for 
improved  absorption  at  boundaries  of  the  simulation  domain. 
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We  show  the  capacity  of  the  RSG  method  in  modeling  wave  scattering  due  to 
smooth  topography,  mesas  and  canyons.  The  sharp  comers  of  mesas  and  canyons  act  as 
strong  scatters  of  seismic  waves. 

3.5.1.  Advantages  of  modeling  wave  propagation  using  a  rotated-staggered 

grid 

The  RSG  finite  difference  has  a  few  remarkable  features  that  are  advantageous  for 
modeling  an  explosion  in  the  vicinity  of  a  high-contrast  scatterer,  such  as  tunnels.  The 
differences  between  the  standard-staggered  grid  (SSG)  and  the  rotated-staggered  grid 
(RSG)  are  shown  in  Figure  12. 

a.  Media  with  high  contrasts  in  elastic  properties 

In  the  standard  staggered  grid  (SSG),  all  model  parameters  are  located  at  different 
positions  within  a  finite  difference  cell  (Figure  12a).  In  RSG,  all  components  of  one 
physical  property,  such  as  stress  or  elastic  moduli,  are  located  at  only  one  location 
(Figure  1 2b).  Therefore,  the  RSG  method  does  not  average  elastic  properties,  but  the  SSG 
method  must  average  the  properties.  When  modeling  wave  propagation  in  media  with 
high  contrasts  in  elastic  properties,  SSG  can  become  inaccurate  and  unstable,  and  the 
RSG  method  can  handle  this  situation  easily. 

b.  Free  surface  conditions 

In  a  traditional  finite  difference  algorithm,  the  free  surface  boundary  condition  has 
to  be  defined  explicitly.  When  simulating  topography  at  the  free  surface,  one  has  to  use  a 
curvilinear  transform  of  the  grid,  which  makes  programming  very  complicated.  In  RSG, 
the  properties  of  air  (vacuum)  and  earth  and  the  free  surface  are  naturally  incorporated.  In 
other  words,  to  realize  a  free  surface,  we  simply  set  the  elastic  moduli  above  the  free 
surface  to  zero  and  the  density  close  to  zero. 

c.  Numerical  stability  and  dispersion 

The  stability  is  increased  by  a  factor  of  yfi  compared  to  SSG  for  the  velocity- 
stress  RSG  method  in  3D  isotropic  media.  The  SSG  has  lower  dispersion  than  the  RSG 
by  a  factor  of  V3  .  The  RSG  allows  the  larger  calculation  time-step  with  the  same  factor. 
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(a) 


(b) 


Figure  12.  Elementary  cells  of  different  staggered  grids  for  a  velocity-stress  finite 
difference  method.  Locations  where  stress,  velocities,  and  elastic  parameters  are 
defined,  (a)  Velocity-stress  finite  difference  technique  using  a  standard  staggered 
grid,  (b)  Velocity-stress  finite  difference  technique  using  the  rotated  staggered  grid. 
(Saenger  and  Bohlen,  2004) 


Therefore,  the  newly  implemented  RSG  finite-difference  program  is  a  very  useful 
tool  to  study  scattered  waves  due  to  near-source  tunnels  and  surface  wave  scattering  due 
to  topography.  Applying  the  RSG  technique  to  a  velocity-stress  formulation  of  the  wave 
equation,  we  can  simulate  the  propagation  of  seismic  waves  in  the  earth  considering 
attenuation. 

3.5.2.  Scattering  due  to  topography 

We  study  the  effects  of  different  topographic  features  on  propagation  and 
scattering  of  seismic  waves  from  explosion  sources  at  shallow  depth.  In  all  the  following 
studies,  we  use  a  pressure  Kelly  wavelet  with  center  frequency  of  10  Hz  as  an  explosion 
source,  as  shown  in  Figure  13.  The  formation  compressional  and  shear  wave  velocities 
and  density  are  3000  m/s,  1700  m/s,  and  2300  kg/m3,  respectively. 
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Figure  13.  A  pressure  Kelly  wavelet  with  center  frequency  of  10  Hz. 
a.  Scattering  due  to  a  hill 

We  first  study  the  effect  of  a  hill  on  wave  propagation.  The  hill  is  300  m  high  and 
2000  m  wide.  A  tunnel  is  placed  100  meters  below  the  free  surface,  as  shown  in  Figure 
14.  An  explosion  source  with  a  center  frequency  of  10  Hz  is  also  placed  at  the  same 
depth.  The  length  and  radius  of  the  cylindrical  tunnel  are  50  m  and  15  m,  respectively. 
The  axis  of  the  tunnel  is  parallel  to  the  y  axis  of  the  Cartesian  coordinates.  The  distance 
between  the  source  and  the  tunnel  axis  is  55  m.  Inside  the  tunnel,  we  assume  an  acoustic 
velocity  of  340  m/s  and  a  compressed  air  density  of  200  kg/m  . 


Figure  14.  Model  for  scattering  due  to  a  hill. 
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Figure  15.  Comparison  of  waveforms  for  a  flat  free  surface  model  and  a  free  surface 
with  a  hill. 

Figure  15  compares  waveforms  for  a  flat  free  surface  model  and  a  free  surface 
with  a  hill.  We  see  that  the  presence  of  the  hill  does  not  affect  the  direct  body  waves,  but 
affects  the  surface  waves. 

b.  Scattering  due  to  a  canyon 

To  clearly  show  the  effects  of  topography  on  wave  propagation,  we  compare  the 
snapshots  of  wave  fields  in  earth  with  a  flat  free  surface  (Figure  16a)  with  those  in  earth 
with  a  canyon  and  a  mesa  (Figure  16b,c). 
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Figure  16.  Snapshot  of  the  divergence  (upper  panel)  and  curl  (lower  panel)  of  the 
wavefield  from  an  explosion  in  earth  with  (a)  a  flat  surface,  (b)  a  canyon,  (c)  a  mesa. 
Note  the  strong  scattering  due  to  the  corners  of  the  features. 
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We  show  that  the  RSG  method  can  effectively  model  the  scattering  due  to  a 
canyon.  The  canyon  model  is  50  m  in  depth  and  width.  An  explosion  source  with  a  center 
frequency  of  10  Hz  is  also  placed  at  100  m  depth  and  200  m  away  from  the  center  of  the 
canyon.  Snapshots  of  the  divergence  and  curl  of  the  wavefield  (Figure  16b)  show  the 
scattering  effects  of  the  sharp  corners  of  the  canyon. 

c.  Scattering  due  to  a  mesa 

The  two-dimensional  mesa  model  is  50  m  high  and  100  wide.  An  explosion 
source  with  a  center  frequency  of  10  Hz  is  also  placed  100  m  below  the  free  surface  and 
200  m  away  from  the  center  of  the  mesa.  Divergence  and  curl  of  the  displacement  field 
separate  the  compressional  and  shear  components  of  the  wavefield.  Figure  16c  shows  the 
snapshots  of  the  divergence  and  curl  of  the  velocity  field.  Complex  patterns  of  the 
scattering  can  be  observed  due  the  existence  of  the  free  surface  and  the  sharp  comers  of 
the  mesa. 


4.  Wavelet-Domain  Inversion  for  Source  Parameters 


4.1  Review  of  Moment  Tensor  Inversion  Method 

Estimation  of  the  seismic  moment  tensor  and  source-time  function  using 
waveform  inversion  has  been  performed  routinely  by  seismologists  to  study  earthquake 
mechanisms  and  source-time  histories.  Many  methods  have  been  developed  since  Gilbert 
and  Dziewonski  (1975)  used  free  oscillation  data  for  their  inversion  (e.g.,  Langston, 
1981;  Kikuchi  and  Kanamori,  1982;  Sipkin,  1982;  Nabelek,  1984).  Different  types  of 
seismic  waveforms,  such  as  long-period  surface  waves  (McCowan,  1976;  Mendiguern, 
1977)  and  low-frequency  body  wave  data  (Stump  and  Johnson,  1977;  Langston,  1981), 
were  inverted  for  source  mechanisms.  However,  there  has  been  limited  success  (e.g., 
Sileny  et  al.,  1992;  Sileny  and  Psencik,  1995;  Schurr  and  Nabelek,  1999)  in  applying 
these  w'aveform  inversion  techniques  to  high-frequency  seismograms,  despite  widespread 
availability  of  broadband  three-component  data.  The  waveform  inversion  methods  being 
used  are  performed  in  either  the  time  or  frequency  domain. 
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In  earthquake  mechanism  waveform  inversions,  the  most  commonly  used  method 
to  reduce  the  number  of  parameters  to  be  determined  is  using  boxcar  functions  (Langston, 
1981)  or  overlapping  triangles  (Nabelek,  1984).  These  parameterization  bases  sometimes 
perform  poorly  in  approximating  the  model  (e.g.,  source-time  function,  moment  tensor 
rate  function,  etc.).  If  these  bases  failed  to  represent  the  model  properly,  they  would  give 
poor  estimates  of  source  parameter.  Therefore,  it  is  preferable  to  choose  a  basis  that  can 
construct  precise  approximations  with  a  linear  combination  of  a  small  number  of  vectors 
selected  inside  the  basis. 

4.2  Wavelet-Domain  Inversion  Method 

In  this  study,  we  will  use  a  wavelet-based  approach  to  formulate  the  waveform 
inversion  problem.  Wavelet  analysis  is  an  increasingly  popular  tool  for  numerical  studies 
in  signal  processing  (Mallat,  1989;  Wang  et  al .,  1995),  biomedical  applications  (Delaney 
and  Bresler,  1995;  Zhu  et  al.,  1997),  and  geophysics  (Deighan  and  Watts,  1997;  Anant 
and  Dowla,  1997;  Wood,  1999;  Kane  and  Herrmann,  2001;  Kane  et  al.,  2002;  Kane, 
2003).  An  extensive  description  of  geophysical  applications  can  be  found  in  Foufoula- 
Georgiou  and  Kumar  (1995),  while  a  theoretical  treatment  of  wavelet  analysis  is  given 
by  Daubechies  (1992),  Strang  and  Nguyen  (1997),  and  Mallat  (1998).  However,  no 
previous  attempt  to  apply  wavelet  analysis  to  earthquake  source  mechanism  inversion  has 
been  undertaken.  In  our  proposed  method,  we  adopt  a  wavelet-based  strategy  to 
parameterize  the  moment  tensor  rate  functions  (MTRFs).  The  MTRFs  allow  the  time- 
dependent  source  mechanism  such  that  each  moment  tensor  component  has  its  own  time 
history  (Dziewonski  and  Gilbert,  1974;  Stump  and  Johnson,  1977;  Ruff  and  Tichelaar, 
1990;  Sileny  et  al.,  1992).  By  choosing  the  “best”  wavelet  as  the  basis,  we  can  construct 
an  adaptive,  problem-dependent  parameterization  for  the  MTRFs,  thereby  achieving 
accurate  approximations  while  significantly  reducing  the  number  of  parameters  that  need 
to  be  estimated  through  inversion.  Additionally,  we  will  perform  inversion  in  the  wavelet 
domain  rather  than  the  time  or  frequency  domain.  This  gives  the  advantage  of  solving 
the  inverse  problem  in  a  multi-resolution  sparse  matrix  representation.  We  can  then  solve 
the  problem  from  coarse  to  fine  levels  out  to  the  limit  of  stability  (Sze  and  Toksoz,  2002). 
At  each  resolution  level,  a  regularized  least-squares  solution  is  obtained  using  the 
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conjugate  gradient  method.  By  significantly  reducing  the  number  of  parameters  that  need 
to  be  estimated,  and  solving  the  inverse  problem  from  coarse  to  fine  levels,  the  wavelet- 
based  method  allows  us  to  obtain  stable  and  geologically  sensible  solutions  more  easily. 
Another  advantage  of  transforming  the  inverse  problem  to  the  wavelet  domain  is  that 
wavelets  are  powerful  tools  for  denoising  data  (Donoho,  1993).  Transforming  data  to  the 
wavelet  domain  tends  to  isolate  signals  into  a  few  large-valued  coefficients,  while  the 
background  noise  tends  to  spread  around  equally  with  less  energy.  Therefore,  we  can 
incorporate  a  nonlinear  wavelet  thresholding  operator  to  remove  the  small  wavelet 
coefficients,  and  thus  the  noise  is  attenuated  with  little  effect  on  the  signals. 

4.3.  Example  of  Wavelet-Domain  Inversion  Applied  to  an  Explosion  and  Double- 
Couple  Source 

We  calculate  synthetic  seismograms  for  a  combination  of  an  explosion  and  a 
double-couple  source.  Using  these  we  invert  for  a  moment  tensor  using  the  wavelet- 
based  method  and  the  conventional  time-domain  method  using  overlapping  triangular 
parameterization.  We  perform  two  inversion  experiments:  the  first  with  noise-free  data, 
and  the  second  with  correlated  white  noise.  Synthetic  data  are  generated  to  simulate  the 
case  where  an  explosive  source  is  accompanied  by  an  earthquake  with  strike-slip 
mechanism.  Figures  17  and  18  show  the  location  of  the  source  and  station  network. 
Assuming  that  the  sources  of  the  explosion  and  the  earthquake  are  at  the  same  place  with 
depth  of  0.5  km  and  have  the  same  source-time  function  with  duration  of  0.5  s.  We 
generate  synthetic  seismograms  with  a  sampling  rate  of  125  Hz  for  both  the  explosive 
and  strike-slip  sources  using  the  reflectivity  method  (Kennett,  1983)  and  then  sum  them. 
We  assume  that  the  duration  of  the  MTRFs  is  one  second  (125  samples).  Results  of  the 
inversion  (Figure  19)  show  that  the  wavelet  method  has  no  difficulty  in  resolving  the 
combined  source  mechanism.  Most  importantly,  we  can  obtain  a  good  of  MTRFs  with 
only  eight  wavelet  coefficients.  In  order  to  compare  our  wavelet-based  method  to  the 
conventional  method,  we  repeat  the  inversion  with  the  same  dataset  in  the  time-domain 
with  eight  overlapping  triangles  for  parameterization.  The  results  are  shown  in  Figure 
20.  It  can  be  seen  that  the  conventional  method  fails  to  recover  a  significant  amount  of 
the  isotropic  components  of  the  MTRFs.  It  also  underestimates  the  magnitude  of  all  the 
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components.  To  achieve  almost  the  same  accuracy,  the  conventional  time  domain  method 
needs  32  overlapping  triangles,  4  times  that  needed  by  wavelet  domain  method. 

In  the  second  experiment,  we  evaluate  the  robustness  of  the  method.  We  add 
correlated  noise  that  is  simulated  by  sampling  an  autoregressive  process  of  order  2.  A 
high  noise  level  is  set  to  be  high,  with  amplitude  comparable  to  the  amplitude  of  the 
signal.  Results  obtained  from  the  wavelet-based  method  and  the  conventional  method  are 
shown  in  Figures  21  and  22,  respectively.  With  the  convenient  implementation  of  the 
wavelet  thresholding  technique,  the  wavelet-based  method  again  outperforms  the 
conventional  method  and  recovers  the  solution  reasonably  well.  Although  the 
conventional  method  is  able  to  recover  part  of  the  isotropic  components  of  the  MTRFs, 
significant  artifacts  show  up  in  the  latter  part  of  the  solution  from  0.5  to  1  sec,  while 
again  underestimating  the  magnitude  of  the  MTRFs. 

The  results  obtained  using  the  wavelet-based  method  show  promise.  We  believe 
this  method  can  be  applied  to  the  understanding  of  the  complexity  of  source  mechanisms 
associated  with  a  nuclear  explosion,  and  has  the  ability  to  resolve  the  volumetric 
components  from  records  that  are  contaminated  by  earthquakes  or  secondary  sources, 
such  as  radiation  from  strong  scatterers. 


29 


100 


90 


80 


70 


E  60 


50 


40- 


30 


20 


20 


40  60 

km 


80 


100 


Figure  17.  Station  positions  for  the  synthetic  seismograms.  All  of  them  are  of 
regional  distances. 
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Figure  18.  Three-component  synthetic  seismogram  of  station  1  generated  by 
reflectivity  method.  It  is  obtained  by  summing  two  seismograms,  one  generated  by 
an  explosive  source  and  the  other  by  a  strike-slip  source  with  a  moment  tensor 
component  of  M23. 
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Figure  19,  Inversion  results  of  the  wavelet-based  method  on  noise-free  data.  Each 
trace  represent  a  component  of  the  MTRFs  (from  top  to  bottom:  the  diagonal 
components  MTU,  MT22,  MT33,  and  the  off-diagonal  components  MT12,  MT13, 
MT23).  The  pink  lines  show  the  original  source-time  functions  as  inputs  to  the 
forward  modeling.  The  black  lines  show  the  MTRFs  recovered  by  the  inversion. 
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Figure  20.  Inversion  results  of  the  conventional  time-domain  method  on  noise-free 
data.  The  pink  lines  show  the  original  source-time  functions  as  inputs  to  the  forward 
modeling.  The  black  lines  show  the  MTRFs  recovered  by  the  inversion. 
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Figure  21.  Inversion  results  of  the  wavelet-based  method  on  synthetic  data 
contaminated  by  correlated  white  noise. 
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Figure  22.  Inversion  results  of  the  conventional  time-domain  method  on  synthetic 
data  contaminated  by  correlated  white  noise. 

4.4  Application  of  Wavelet-Domain  Inversion  for  Explosion  near  a  Tunnel 

We  use  synthetic  seismograms  to  test  the  performance  of  the  wavelet-domain 
moment  tensor  inversion  method  on  seismograms  generated  by  the  finite  difference 
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approach.  The  schematic  diagram  of  the  explosion,  tunnel  and  recording  plane  is  shown 
in  Figure  5.  This  will  enable  us  to  evaluate  how  well  our  method  can  resolve  the  isotropic 
moment  tensor  components  in  the  presence  of  stronger  scatter. 

Figures  23a  and  23b  show  the  results  of  the  wavelet-domain  moment  tensor 
inversion  and  the  fitting  of  the  waveforms  for  the  explosion  with  a  tunnel,  respectively. 
The  wavelet-domain  inversion  method  was  able  to  retrieve  the  explosive  source 
mechanism,  indicated  by  the  large  isotropic  components  and  zero  deviatoric  components 
before  0.05  s.  The  effect  of  scattering  from  the  tunnel  only  introduced  an  anomalous,  late 
shear  event  at  around  0. 15  s  (Figure  23a)  and  did  not  affect  our  ability  to  resolve  the 
original  explosive  mechanism.  We  obtained  good  fittings  of  the  synthetic  waveforms  as 
well. 
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Figure  23a:  The  results  of  wavelet-domain  moment  tensor  inversion  for  an  explosion 
with  tunnel.  Some  anomalous  components  of  Mil,  M22  and  M13  showed  up  at 
a  later  time  (-0.15  s)  but  the  inversion  was  still  able  to  retrieve  large  isotropic 
components  (Mil,  M22,  M33)  and  zero  deviatoric  components  (M12,  M13, 
M23)  in  the  beginning. 
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Figure  23b:  Fitting  of  waveforms  (blue)  by  the  wavelet-domain  inversion  to  the 
synthetic  data  (black)  of  an  explosion  with  tunnel.  Data  from  eight  stations 
(circles)  surrounding  the  explosion  (cross)  were  used  for  the  inversion.  At  each 
station  traces  are  (top  to  bottom)  vertical,  tangential  and  radial  components. 


5.  Source  and  Scatter  Imaging  Using  Time  Reversed  Acoustics 


5.1  Introduction  of  Time  Reversed  Acoustics 

According  to  the  TRA  concept,  acoustic  waves  recorded  at  several  stations  when 
time-reversed  and  put  back  into  the  medium,  propagate  and  focus  at  the  original  source 
(Fink,  1993,  1999,  2001;  Song  and  Kupperman,  1999;  Derode  et  al.,  2000;  Lu  and 
Toksoz,  2004).  The  concept  is  illustrated  schematically  in  Figures  24a, b.  To 
demonstrate  the  applicability  to  seismic  source  characterization,  we  conducted  two 
numerical  experiments.  In  each  experiment,  seismic  waves  generated  by  a  source  in 
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heterogeneous  medium  were  calculated  using  the  finite  difference  code.  Then  these 
seismograms  were  time  reversed  and  “pumped”  back  into  the  medium. 

5.2  Explosion  in  a  layered  Homogeneous  Medium 

Figure  25  shows  the  geometry  of  the  first  experiment  where  an  explosion  source 
is  placed  in  a  layered  medium  and  recorded  by  a  circular  array  of  receivers.  Note  that  the 
source  is  off-centered.  Figure  26  shows,  through  a  series  of  snapshots,  the  convergence 
of  the  waves  to  the  source. 

5.3  Explosion  near  a  Tunnel  Buried  in  a  Homogeneous  Medium 

The  next  example  is  for  an  explosion  that  is  located  near  a  tunnel.  A  subsurface 
explosion  is  modeled  by  the  elastic  finite  difference  code  and  the  seismograms  generated 
are  recorded  at  multiple  surface  stations.  The  TRA  methodology  is  applied  to  these 
recorded  seismograms  by  first  reversing  them  in  time  and  then  applying  them  as  sources 
at  the  top  of  the  same  3D  elastic  modeling  algorithm.  Figure  27  shows  snapshots  of  the 
back  propagated  wave  field  at  four  different  times.  The  left  column  shows  the  horizontal 
component  of  motion,  the  middle  column  shows  the  divergence  of  the  wave  field  or  P 
wave  only  energy,  and  the  right  column  shows  the  curl  of  the  wave  field  or  the  shear 
wave  energy  only.  The  star  indicates  the  location  of  the  source  and  the  circle  shows  the 
outline  of  the  tunnel. 

The  record  length  of  the  data  used  for  back  propagation  was  0.4595  s.  The  top 
row  of  snapshots  represent  0.2995  s  into  the  back  propagating  model  which  corresponds 
to  0.160  s  in  the  original  forward  model  time.  In  a  similar  fashion,  the  bottom  row  of 
snapshots  represent  0.4595  s  into  the  back  propagating  model  or  0.0  s  in  the  original 
model  time.  Clearly  seen  in  the  bottom  row  at  the  original  model  time  is  the  convergence 
of  the  P  wave  energy  at  the  source  position  in  the  left  and  middle  figures.  The  shear 
wave  energy  converges  to  the  tunnel,  the  source  of  the  P  to  shear  wave  conversion,  at  a 
back  propagation  time  of  0.4405  s. 
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Figure  24a.  Schematic  showing  energy  from  a  seismic  source 
propagating  through  a  medium  with  many  scatterers,  and  being 
recorded  at  the  stations  denoted  by  red  triangles.  The  yellow  traces 
represent  the  recorded  seismograms.  The  recorded  traces  are 
created  by  a  convolution  of  the  source  function,  s(t),  with  the 
appropriate  transfer  function,  gj(t). 
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Figure  24b.  Schematic  showing  the  TRA  process  on  the  data  recorded  as 
in  Figure  4a.  First  the  recorded  seismograms  (top  yellow  traces)  are 
reversed  in  time  (bottom  yellow  traces)  and  pumped  back  into  the 
medium  at  the  station  locations.  The  energy  reverses  its  path 
through  the  medium  and  converges  upon  the  original  source 
position.  Due  to  reciprocity,  the  Greens  function,  gj(t)  for  each 
receiver  to  source  path  is  identical  to  the  source  to  receiver  Greens 
function. 
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Figure  25.  Simple  2D  layered  earth  model  with  seismic  source  located  at  the 
yellow  star,  and  a  ring  of  receivers,  each  indicated  by  a  red  triangle. 


pressure 


Figure  26.  Progressive  snap-shots  of  the  wave  field  propagated  backward  in  time,  as 
originally  recorded  by  the  receivers  in  Figure  25.  The  fourth  column 
corresponds  to  zero  time  in  the  original,  forward  model. 
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Figure  27.  Snapshots  of  back  propagated  records  for  an  explosion  (star)  located  next 
to  an  air  filled  tunnel  (circle).  The  left  column  shows  the  vertical  component  of 
motion,  the  middle  column  shows  the  divergence  of  the  field  (P  wave  energy  only), 
and  the  right  column  shows  the  curl  of  the  field  (shear  energy  only).  The  bottom 
row  represents  zero  time  in  the  original  forward  model. 


6.  Conclusions  and  Recommendations 

1 .  A  tunnel  near  an  explosion  acts  as  a  strong  scatterer.  P  to  S  scattering  is  much  stronger 
than  P  to  P  scattering.  Some  energy  is  scattered  into  SH  waves. 

2.  Surface  topographic  features  act  as  scatterers.  Overall,  scattering  from  topography  is 
small  compared  to  that  of  tunnels.  The  smooth  topography  effect  on  scattering  is 
relatively  small.  Sharp  topographic  features  cause  strong  scattering  of  the  surface  waves. 

3.  For  explosion  identification,  wavelet-based  moment  tensor  inversion  recovers  the 
relative  strengths  of  the  explosion  and  the  scatterer  source. 
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4.  TRA  methodology  has  the  potential  of  determining  seismic  source  properties  with 
good  resolution. 
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